{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from classy import Class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_out = [5e-5, 5e-4, 5e-3]\n",
    "models = ['PPF1','PPF2','FLD1','FLD1S']\n",
    "w0 = {'PPF1':-0.7,'PPF2':-1.15,'FLD1':-0.7,'FLD1S':-0.7}\n",
    "wa = {'PPF1':0.,'PPF2':0.5,'FLD1':0.,'FLD1S':0.}\n",
    "omega_cdm = {'PPF1':0.104976,'PPF2':0.120376,'FLD1':0.104976,'FLD1S':0.104976}\n",
    "omega_b = 0.022\n",
    "##Omega_cdm = {'PPF1':0.26,'PPF2':0.21,'FLD1':0.26,'FLD1S':0.26}\n",
    "##Omega_b = 0.05\n",
    "h = {'PPF1':0.64,'PPF2':0.74,'FLD1':0.64,'FLD1S':0.64}\n",
    "cosmo = {}\n",
    "\n",
    "for M in models:\n",
    "    use_ppf = 'yes'\n",
    "    gauge = 'Newtonian'\n",
    "    if 'FLD' in M:\n",
    "        use_ppf = 'no'\n",
    "    if 'S' in M:\n",
    "        gauge = 'Synchronous'\n",
    "        \n",
    "    cosmo[M] = Class()\n",
    "    \n",
    "    cosmo[M].set({'output':'tCl mPk dTk vTk','k_output_values':str(k_out).strip('[]'),\n",
    "                  'h':h[M],\n",
    "                  'omega_b':omega_b,'omega_cdm':omega_cdm[M],\n",
    "                  ##'Omega_b':Omega_b,'omega_cdm':Omega_cdm[M],\n",
    "                  'cs2_fld':1.,\n",
    "          'w0_fld':w0[M],'wa_fld':wa[M],'Omega_Lambda':0.,'gauge':gauge,\n",
    "                 'use_ppf':use_ppf})\n",
    "    cosmo[M].compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colours = ['r','k','g','m']\n",
    "for i,M in enumerate(models):\n",
    "    cl = cosmo[M].raw_cl()\n",
    "    l = cl['ell']\n",
    "    \n",
    "    plt.loglog(l,cl['tt']*l*(l+1)/(2.*np.pi),label=M,color=colours[i])\n",
    "    \n",
    "plt.legend(loc='upper left')\n",
    "plt.xlim([2,300])\n",
    "plt.ylim([6e-11,1e-9])\n",
    "plt.xlabel(r'$\\ell$')\n",
    "plt.ylabel(r'$[\\ell(\\ell+1)/2\\pi]  C_\\ell^\\mathrm{TT}$')\n",
    "\n",
    "plt.savefig('check_PPF_clTT.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for M in ['PPF1','FLD1']:\n",
    "    csm = cosmo[M]\n",
    "    pt = csm.get_perturbations()\n",
    "    pts = pt['scalar']\n",
    "    for i,k in enumerate(k_out):\n",
    "        ptk = pts[i]\n",
    "        a = ptk['a']\n",
    "        phi = ptk['phi']\n",
    "        psi = ptk['psi']\n",
    "        if 'FLD' in M:\n",
    "            ls = ':'\n",
    "            lw=5\n",
    "        else:\n",
    "            ls = '-'\n",
    "            lw=1\n",
    "        plt.semilogx(a,0.5*(phi+psi),label=M+' '+'$k='+str(k)+'Mpc^{-1}$',ls=ls,lw=lw)\n",
    "        \n",
    "plt.legend(loc='lower left')\n",
    "plt.xlim([1e-2,1])\n",
    "plt.ylim([0.3,0.63])\n",
    "plt.xlabel(r'$a/a_0$')\n",
    "plt.ylabel(r'$\\frac{1}{2} ~(\\Phi+\\Psi)$')\n",
    "\n",
    "plt.savefig('check_PPF_metric.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#kminclosed = sqrt(-8*Omega_k)*(70/3e5) Mpc^(-1)\n",
    "\n",
    "k_out = [1e-3] #[1e-4, 1e-3, 1e-2]\n",
    "#models = ['PPF1','PPF2','FLD1']\n",
    "models = ['PPF1','FLD1']\n",
    "w0 = {'PPF1':-0.7,'PPF2':-1.15,'FLD1':-0.7,'FLD1S':-0.7}\n",
    "wa = {'PPF1':0.,'PPF2':0.5,'FLD1':0.,'FLD1S':0.}\n",
    "omega_cdm = {'PPF1':0.104976,'PPF2':0.120376,'FLD1':0.104976,'FLD1S':0.104976}\n",
    "omega_b = 0.022\n",
    "##Omega_cdm = {'PPF1':0.26,'PPF2':0.21,'FLD1':0.26,'FLD1S':0.26}\n",
    "##Omega_b = 0.05\n",
    "h = {'PPF1':0.64,'PPF2':0.74,'FLD1':0.64}\n",
    "\n",
    "fig, axes = plt.subplots(1,2,figsize=(16,5))\n",
    "for Omega_K in [-0.1, 0.0, 0.15]:\n",
    "    for gauge in ['Synchronous','Newtonian']:\n",
    "        cosmo = {}\n",
    "        for M in models:\n",
    "            use_ppf = 'yes'\n",
    "            if 'FLD' in M:\n",
    "                use_ppf = 'no'\n",
    "        \n",
    "            cosmo[M] = Class()\n",
    "    \n",
    "            cosmo[M].set({'output':'tCl mPk dTk vTk','k_output_values':str(k_out).strip('[]'),\n",
    "                  'h':h[M],\n",
    "                  'omega_b':omega_b,'omega_cdm':omega_cdm[M],'Omega_k':Omega_K,\n",
    "                  ##'Omega_b':Omega_b,'omega_cdm':Omega_cdm[M],\n",
    "                  'cs2_fld':1.,\n",
    "          'w0_fld':w0[M],'wa_fld':wa[M],'Omega_Lambda':0.,'gauge':gauge,\n",
    "                 'use_ppf':use_ppf,'hyper_sampling_curved_low_nu':10.0})\n",
    "            cosmo[M].compute()\n",
    "            \n",
    "        label = r'$\\Omega_k='+str(Omega_K)+'$, '+gauge[0]\n",
    "        clfld = cosmo['FLD1'].raw_cl()\n",
    "        clppf = cosmo['PPF1'].raw_cl()\n",
    "        \n",
    "        axes[0].semilogx(clfld['ell'][2:],clppf['tt'][2:]/clfld['tt'][2:],label=label)\n",
    "        \n",
    "        ptfld = cosmo['FLD1'].get_perturbations()['scalar']\n",
    "        ptppf = cosmo['PPF1'].get_perturbations()['scalar']\n",
    "        for i,k in enumerate(k_out):\n",
    "            ptkfld = ptfld[i]\n",
    "            a = ptkfld['a']\n",
    "            phi_plus_phi_fld = ptkfld['phi']+ptkfld['psi']\n",
    "            ptkppf = ptppf[i]\n",
    "            phi_plus_phi_ppf = ptkppf['phi']+ptkppf['psi']\n",
    "            axes[1].semilogx(ptkppf['a'],phi_plus_phi_ppf,label=label+'_ppf')\n",
    "            axes[1].semilogx(ptkfld['a'],phi_plus_phi_fld,label=label+'_fld')\n",
    "            print (len(ptkppf['a']),len(ptkfld['a']))\n",
    "            \n",
    "axes[0].legend(loc='lower left',ncol=2)\n",
    "axes[0].set_xlim([2,300])\n",
    "axes[0].set_ylim([0.98,1.02])\n",
    "axes[0].set_xlabel(r'$\\ell$')\n",
    "axes[0].set_ylabel(r'$C_\\ell^\\mathrm{FLD1}/C_\\ell^\\mathrm{PPF1}$')\n",
    "\n",
    "axes[1].legend(loc='lower left',ncol=2)\n",
    "axes[1].set_xlim([1e-2,1])\n",
    "axes[1].set_xlabel(r'$a/a_0$')\n",
    "axes[1].set_ylabel(r'$(\\Phi+\\Psi)$')\n",
    "\n",
    "fig.savefig('check_PPF_Omegak.pdf')        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colours = ['r','k','g','m']\n",
    "\n",
    "k_out = [1e-1] #[1e-4, 1e-3, 1e-2]\n",
    "#models = ['PPF1','PPF2','FLD1']\n",
    "models = ['PPF1','FLD1']\n",
    "w0 = {'PPF1':-0.7,'PPF2':-1.15,'FLD1':-0.7,'FLD1S':-0.7}\n",
    "wa = {'PPF1':0.,'PPF2':0.5,'FLD1':0.,'FLD1S':0.}\n",
    "omega_cdm = {'PPF1':0.104976,'PPF2':0.120376,'FLD1':0.104976,'FLD1S':0.104976}\n",
    "omega_b = 0.022\n",
    "##Omega_cdm = {'PPF1':0.26,'PPF2':0.21,'FLD1':0.26,'FLD1S':0.26}\n",
    "##Omega_b = 0.05\n",
    "h = {'PPF1':0.64,'PPF2':0.74,'FLD1':0.64}\n",
    "\n",
    "fig, axes = plt.subplots(1,2,figsize=(18,8))\n",
    "\n",
    "for Omega_K in [-0.1, 0.0, 0.15]:\n",
    "    for ppfgauge in ['Synchronous','Newtonian']:\n",
    "        cosmo = {}\n",
    "        for M in models:\n",
    "            use_ppf = 'yes'\n",
    "            gauge = ppfgauge\n",
    "            if 'FLD' in M:\n",
    "                use_ppf = 'no'\n",
    "        \n",
    "            cosmo[M] = Class()\n",
    "    \n",
    "            cosmo[M].set({'output':'tCl mPk dTk vTk','k_output_values':str(k_out).strip('[]'),\n",
    "                  'h':h[M],\n",
    "                  'omega_b':omega_b,'omega_cdm':omega_cdm[M],'Omega_k':Omega_K,\n",
    "                  ##'Omega_b':Omega_b,'omega_cdm':Omega_cdm[M],\n",
    "                  'cs2_fld':1.,\n",
    "          'w0_fld':w0[M],'wa_fld':wa[M],'Omega_Lambda':0.,'gauge':gauge,\n",
    "                 'use_ppf':use_ppf,'hyper_sampling_curved_low_nu':6.1})\n",
    "            cosmo[M].compute()\n",
    "            \n",
    "        #fig, axes = plt.subplots(1,2,figsize=(16,5))\n",
    "        for j,M in enumerate(models):\n",
    "            cl = cosmo[M].raw_cl()\n",
    "            l = cl['ell']\n",
    "            label = M+r'$\\Omega_k='+str(Omega_K)+'$, '+gauge[0]\n",
    "            axes[0].loglog(l,cl['tt']*l*(l+1)/(2.*np.pi),label=label,color=colours[j])\n",
    "        \n",
    "            csm = cosmo[M]\n",
    "            pt = csm.get_perturbations()\n",
    "            pts = pt['scalar']\n",
    "            for i,k in enumerate(k_out):\n",
    "                ptk = pts[i]\n",
    "                a = ptk['a']\n",
    "                phi = ptk['phi']\n",
    "                psi = ptk['psi']\n",
    "                if 'FLD' in M:\n",
    "                    ls = ':'\n",
    "                    lw=5\n",
    "                else:\n",
    "                    ls = '-'\n",
    "                    lw=1\n",
    "                axes[1].semilogx(a,0.5*abs(phi+psi),label=label+' '+'$k='+str(k)+'Mpc^{-1}$',ls=ls,lw=lw)\n",
    "\n",
    "axes[0].legend(loc='upper left')\n",
    "axes[0].set_xlim([2,300])\n",
    "axes[0].set_ylim([6e-11,1e-9])\n",
    "axes[0].set_xlabel(r'$\\ell$')\n",
    "axes[0].set_ylabel(r'$[\\ell(\\ell+1)/2\\pi]  C_\\ell^\\mathrm{TT}$')\n",
    "\n",
    "axes[1].legend(loc='upper right')\n",
    "#axes[1].set_xlim([1e-2,1])\n",
    "#axes[1].set_ylim([0.3,0.63])\n",
    "axes[1].set_xlabel(r'$a/a_0$')\n",
    "axes[1].set_ylabel(r'$\\frac{1}{2}~(\\Phi+\\Psi)$')\n",
    "\n",
    "fig.savefig('check_PPF_Omegak2.pdf')               "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (0.31*0.64**2-0.022)\n",
    "print (0.26*0.74**2-0.022)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
